Traffic jams, gliders, and bands in the quest for collective motion 
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We study a simple swarming model on a two-dimensional lattice where the self-propelled particles 
exhibit a tendency to align ferromagnetically. Volume exclusion effects are present: particles can 
only hop to a neighboring node if the node is empty. Ifere we show that such effects lead to 
a surprisingly rich variety of self-organized spatial patterns. As particles exhibit an increasingly 
higher tendency to align to neighbors, they first self-segregate into disordered particle aggregates. 
Aggregates turn into traffic jams. Traffic jams evolve toward gliders, triangular high density regions 
that migrate in a well-defined direction. Maximum order is achieved by the formation of elongated 
high density regions - bands - that transverse the entire system. Numerical evidence suggests that 
below the percolation density the phase transition associated to orientational order is of first-order, 
while at full occupancy it is of second-order. 
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Self-propelled particle (SPP) systems are found at all 
scales in nature. Examples in biology range from human 
crowds [3l and animal groups down to insects 

bacteria [6], and even to the microcellular scale with, e.g., 
the collective motion of microtubules driven by molecu- 
lar motors [7]. SPP systems are not restricted to living 
systems. There are examples in non-living matter, as for 
instance in driven granular media [sl-fiot. Interestingly, 
the statistical properties of the large-scale self-organized 
patterns emerging in SPP systems depend only on a few 
microscopic details: the symmetry associated to the self- 
propulsion mechanism of the particles, which can be ei- 



ther polar llHlSj or apolar the symmetry of the 



velocity alignment mechanism, which can be either fer- 

isl 14 1, and very impor- 



11, 12 or nematic 



romagnetic 

tantly, the presence or absence of volume exclusion ef- 
fects, as we are going to discuss here. In addition, the 
nature of the supporting space where the particles move 
plays also a crucial role: the dimension of the space, and 
whether this space is continuous 11-14| or discrete [lil - 

a. 



In this letter, we focus on polar SPPs moving on a 
two-dimensional lattice that align their orientation, re- 
spectively their moving direction, via a local ferromag- 
netic alignment mechanism. We explicitly model vol- 
ume exclusion effects: nodes can be occupied at most 
by one particle. We show that such effects introduce a 
coupling between particle speed, local density and local 
alignment that lead to a surprisingly rich variety of self- 
organized spatial patterns unseen in previous swarming 
models. As particles exhibit an increasingly higher ten- 
dency to align to neighbors, the system passes through 
three distinct phases. For weak alignment strength, the 
system exhibits orientational disorder, while particles 
self-segregate. Fig. [IJa). Within this initial phase, there 
is a transition from a spatially homogeneous to an ag- 



gregate phase. The onset of orientational (polar) order 
marks the beginning of the second phase which is charac- 
terized by the emergence of locally ordered, high density 
regions: traffic jams, Fig. [TJb). As the tendency to align 
is enhanced, traffic jams evolve toward triangular high 
density aggregates that migrate in a well-defined direc- 
tion. We refer to these dynamical traffic jams as gliders. 
Fig. [He). The third phase emerges when the particles 
self-organize into highly ordered, elongated, high density 
regions: bands. Fig. [Ijd). In contrast to the traveling 
bands observed in off-lattice SPP models with ferromag- 
netic alignment jl^l , these bands are formed by particles 
aligned to the long axis of the band and are rather static. 

We find evidence that the phase transition to ori- 
entational order is discontinuous below the percolation 
threshold. When the lattice is fully occupied, the sys- 
tem reduces to the classical planar Potts model and the 
phase transition to orientational order is undoubtedly of 
second-order |2J, |25|. Previous lattice swarming mod- 



els were found to exhibit a continuous phase transition 
from a homogeneous to a condensed phase in ID (l7l.[l8j. 
while in 2D, both, first and second-order transitions to 
orientational order have been claimed. Bussemaker et al. 
reported a second-order transition in a cellular automa- 
ton model with 4 moving directions 15[, while Csahok 
and Vicsek found, for a lattice-gas model with 6 mov- 
ing directions, a weakly-first-order transition to collec- 
tive migration [l^ . Here we show that all these phenom- 
ena occur in a minimal 2D lattice swarming model, but 
where in contrast to previous models the system dynam- 
ics is dominated by volume exclusion effects that lead to 
a completely novel spatial self-organization of particles. 

Model. Particles move on a two-dimensional lattice with 
periodic boundary conditions, and have four possible ori- 
entations: up, down, left, and right, whose associated 
vectors are vi, V2, V3, and V4, respectively. The state of 
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FIG. 1: (Color online) Example of self-organized spatial patterns: (a) disordered aggregates, (b) traffic jams, (c) gliders, and 
(d) bands. Particle orientation is indicated by the orientation of the small triangles, and is also color-coded: right (red), left 
(black), down (green), and up (blue). Parameter values are indicated in Fig. [2l 



a particle is given by its position on the lattice and its 
orientation. As it will become clear below, the orienta- 
tion of a particle fully determines its moving direction. 
The particles are able to perform two actions: i) they 
can change their orientation, and ii) they can migrate in 
the direction given by their orientation. These actions 
have associated transition rates which specify the aver- 
age number of events per time unit. Let us start out with 
the reorientation transition rate Tr that a particle at x 
with orientation v turns its orientation into direction w: 

r;^((x,v)^(x,w)) = exp(5 (^Iv(y)», (1) 

y 6A(x) 

where the sum runs over the nearest lattice neighbors of 
X, represented by A(x). The vector V(y) returns the 
orientation of the particle placed in node y, if there is 
any, and the null vector otherwise. The symbol (.|.) in- 
dicates the inner product between two vectors, while g 
is a parameter which controls the alignment sensitivity. 
For positive g, eq. ([T]) defines a stochastic ferromagnetic 
alignment mechanism. As result of this alignment, first 
nearest neighbor particles tend to be aligned. 

The migration rate is defined in the following way: 

Tm ((x, v) (y = X + V, v)) 

_ ( vq if node y is empty (2) 
[0 if node y is occupied. 

Thus, a particle at position x and pointing in direction 
V migrates to the neighboring node y = x + v with a 
transition rate vq as long as the node at y is empty. If 
the y-node is occupied, the particle will not jump. The 
only action that is allowed to the particle in this situation 
is to change its orientation. 

At this point, it is important to understand how this 
continuous-time process is simulated. Let us assume that 
at time to the system is in a given state. Then, we com- 
pute the time at which the next event will take place in 
the system, i.e., we calculate ti. Now we have to de- 
cide which is the event that will take place at time ti. 
We choose at random one out of all the possible events 
that could take place, but we weight each of these events 
according to their associated transition rate, eqs. ([T|) 




FIG. 2: (Color online) Mean orientation (m) vs. the align- 
ment sensitivity g for systems with density d = 0.3 and mi- 
gration rate vo = 100. Simulations were carried out for 2 10^ 
time steps. The different curves correspond to different sys- 
tem sizes L. The boundary between phase I and II, g* , and 
the between phase II and III, gt are indicated by the verti- 
cal dashed lines, (a) to (d) refer to the simulation snapshots 
shown in Fig. [T] 



and This procedure is an adaption of the classi- 

cal Gillespie algorithm [19|] to interacting particle sys- 
tems |20|. 

Results. We start by fixing the migration rate vq 
and particle density d, and use g as control parame- 
ter. The degree of orientational order in the system 
at time t is characterized by the (global) orientation 
m{t) = {l/N)\J2^Y{x)\, where is the total number 
of particles in the system, the sum runs over all lattice 
sites, and V(x) is defined as above. Fig. [5] shows the 
behavior of the mean orientation (m) as function of the 
alignment sensitivity g, with (...) a temporal average 
taken once the system reaches the steady-state after a 
short transient. The system exhibits a phase transition 
to orientational order above a critical g* for all densities, 
as long as vq > 0. Here we focus on < d < dp, where dp 
refers to the (site) percolation threshold in a 2D (square) 
lattice, dp ^ 0.59. The system exhibits three phases with 
g, labeled I, II, and III by increasing alignment sensitiv- 
ity g. Phase I corresponds to g < g* and is characterized 
by exhibiting no macroscopic orientational order. Fig. 
[3] shows that within phase I there is a dynamic phase 
transition from a spatially homogeneous to an aggregate 
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FIG. 3: (Color online) Phase transition to aggregation (phase 
I), (a) to (d) correspond to simulation snapshots whose value 
of g is indicated in (f) (other parameters as in Fig. [2]). The 
average cluster size (A;) normalized by the total number of 
particles, d.L^, as function of g is shown in (e) for various 
system sizes, and in (f) for various densities. 



phase as g is increased. The degree of aggregation is 
characterized by the average cluster size (/c), which is 
computed as the temporal average of (fc) — kp{k,t), 
with p{k,t) = knk(t)/N, where nk{t) is the number of 
clusters of mass k at time t and N = dL'^ is the number 
of particles in the system. The figure shows that there 
exists a critical value ga above which a phase transition 
to aggregation occurs, with ga < g* ■ g approaches 
g*, more than 85% of the particles in the system form 
a large aggregate. The transition point between phase I 
and II is at g*, where the curves {m){g,L) for different 
system sizes L meet. Interestingly, g* seems to be inde- 
pendent of the density d, as confirmed with simulations 
for various system sizes with density d = 0.2, 0.3, and 
0.4 (data not shown). Moreover, g* coincides with the 
critical point for the full occupancy problem, i.e., d = 1. 
Fig. Ufa) shows time series of m(t) for values of g close 
to g*. The order parameter m(t) exhibits fluctuations 
between high and low values. Low m-values correspond 
to the appearance of round traffic jams, while high values 
correspond to elongated traffic jams where two directions 
dominate over the other two. Traffic jams results from 
the jamming of four particle clusters attempting to move 
to the left, right, upward, and downward, respectively. 
Fluctuations are due to the competition between these 
four clusters. Fig.lU^b) shows that the distributions p(m) 



FIG. 4: (Color online) Phase transition to orientational or- 
der. Time series m{t), histogram, and typical spatial configu- 
rations for phase II, g* < g < gt- Letters in (c) indicate sim- 
ulation snapshots shown in (d)-(f). Parameters as in Fig. [21 



obtained from the m{t) time series for values of g close 
to g* do not exhibit a Gaussian shape as expected for 
a second order transition, but rather a bimodal distribu- 
tion (the arrow indicates the second peak) as expected for 
first order transitions. The coexistence of several particle 
configurations (or "phases") exhibiting different degrees 
of ordering, i.e., values of to, is evident for values of g 
deep into phase II. Fig. \^c) shows that m{t) jumps be- 
tween well-defined values corresponding to different spa- 
tial patterns. Values of m{t) ~ 1/a/2 are associated to 
gliders, while higher values of m{t) correspond to bands. 
Lower values of m{t) are due to the presence of traffic 
jams. Gliders are dynamical traffic jams moving back- 
wards with respect to their mean average orientation. 
Their presence affects the temporal evolution of the cen- 
ter of mass of the system, Xcm(t), which exhibits bal- 
listic motion whenever there is a glider, while otherwise 
is Brownian. As result of this, the average speed of the 
center of mass (Vcm) peaks at values of g where gliders 
are more stable [211, Fig. [S] Gliders are remarkably 
different from traffic jams observed in 2D traffic mod- 



els [24I, arguably due to the presence of the alignment 



mechanism, Eq. ([T]). How frequently gliders appear and 
for how long they survive, depends on the value of g and 
L. For example, for g = 1.4, m{t) displays excursions 
from low to high values that reflect the fact the system 
alternates between traffic jams and gliders. For an illus- 
tration of this dynamics, see 23|. For 5 > g?,, i.e. in 
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FIG. 5: (Color online) (a) average speed of the center of mass 
(Vcm) as function of g. (b) Trajectories of the center of mass 
for two different values of g. 



phase III, the only stable configurations is a band, 
contrast to g*, is highly dependent on L. 



In 



Discussion. In the limit of full occupancy, d — \, par- 
ticles are frozen in their positions and the only action 
allowed to them is reorientation. The system defined by 
Eqs. ([T]) and ^ becomes an equilibrium system in this 
limit, whose order is again characterized by m{t). Since, 
by definition, there are only four possible orientations, 
we can safely claim that the model reduces to a 4-state 
planar Potts model 



24| 



^ It has been shown that this 

model can be reduced to the standard 2-Potts model "25*] . 
In two-dimensions, the standard g-Potts model exhibits 
a continuous transition for q < 4, and a discontinuous 
one for q > 4 f23|. Thus, our model exhibits a second- 
order transition for d — 1, as confirmed via simulations 
(data not shown). For vq > and d < 1, we are in a 
pure non-equilibrium scenario. The migration rule, Eq. 
([5]), breaks detailed balance and prevents us from writing 
down a free energy. Nevertheless, it is worth to compare 
our system with its equilibrium counterpart, the diluted 
Potts model with annealed vacancies. If we represent 
the absence of particles in a lattice position with an ex- 
tra vector direction, we end up with a 5-Potts system 
with full occupancy, instead of a diluted 4-Potts system. 
In consequence, according to what was said above, the 
transition would be first-order. The argument applies to 
the standard Potts model and assumes vacancies are in 
thermal equilibrium, which is not true for Eq. Nev- 
ertheless, it helps to realize that a discontinuous dynamic 
phase transition is quite possible in our non-equilibrium 
system. 

In summary, we have shown through a minimal model 
that volume exclusion effects, when they are allowed to 
stop particle motion, can lead to a surprisingly rich va- 
riety of self-organized patterns. Such effects introduce a 
coupling between local density, local orientation and par- 
ticle speed that strongly affects the large-scale behavior of 
the system, with the jamming of particles playing a dom- 
inant role. This coupling is present in many real systems 
as in gliding bacteria, animal groups, etc. Certainly, sev- 
eral features of the self-organized patterns described here 
depend on the discrete nature of the model. Nevertheless, 



we expect similar phenomena to emerge in off-lattice, 
continuum symmetry systems. For instance, static traf- 
fic jams are probably a robust property of all systems 
where stagnation can occur. Here we have also learned 
that the jamming of self-propelled particles can lead to 
unexpected self-organized structures in two-dimensions 
like dynamical trafhc jams, e.g., gliders. The presence of 
an alignment mechanism induces (local) orientational or- 
der, and provided particles are oriented, density waves of 
stagnated particles should emerge. The results reported 
here are a first step toward a deep understanding of the 
possible phenomena that such a coupling may induce. 
We thank A. Greven and C.F. Lee for useful comments. 



[1] 
[2] 
[3] 

[4; 

[5 
[6 
[7 

[s: 

[9: 
[lo: 

[11 

[12 
[13 

[14 

[15 

[16 
[17 

[18 

[19 
[20 

[21 

[22; 

[23; 

[24' 
[25 



D. Helbing, I. Parkas, and T. Vicsek, Nature (London) 
407, 487 (2000). 

A. Cavagna et al, Proc. Natl. Acad. Sci. 107, 11865 
(2010). 

K. Bhattacharya and T. Vicsek, New J. Phys. 12, 093019 
(2010). 

J. Buhl et al., Science 312, 1402 (2006). 

P. Romanczuk, I.D. Couzin, and L. Schimansky-Geier, 

Phy. Rev. Lett. 102, 010602 (2009). 

H.P. Zhang et al., Proc. Natl. Acad. Sci. 107, 13626 

(2010). 

V. Schaller et al., Nature 467, 73 (2010). 

V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 

105 (2007). 

A. KudroUi et al., Phys. Rev. Lett. 100, 058001 (2008); 
A. KudroUi, Phys. Rev. Lett. 104, 088001 (2010). 
J. Deseigne, O. Dauchot, and H. Chate, Phys. Rev. Lett. 
105, 098001 (2010). 

T. Vicsek et al., Phys. Rev. Lett. 75, 1226 (1995). 

G. Gregoire and H. Chate, Phys. Rev. Lett. 92, 025702 
(2004). 

F. Peruani, A. Deutsch, and M. Bar, Eur. Phys. J-Spec. 
Top. ,157, 111 (2008); F. GineUi et al., Phys. Rev. Lett. 
104, 184502 (2010). 

H. Chate, F. Ginelli and R. Montague, Phys. Rev. Lett. 
96, 180602 (2006). 

H.J. Bussemaker, A. Deutsch, and E. Geigant, Phys. Rev. 
Lett. 78, 5018 (1997). 

Z. Csahok and T. Vicsek, Phys. Rev. E 52, 5297 (1995). 
O.J. O'Loan and M.R. Evans, J. Phys. A: Math. Gen. 
32, 99 (1999). 

J.R. Raymond and M.R. Evans, Phys. Rev. E 73, 036112 
(2006). 

D.T. Gillespie, J. Phys. Chem. 81, 2340 (1977). 
T. Klauss and A. Voss-Boehme, in Mathematical Model- 
ing of Biological Systems, vol. II, 353 (2008). 
(Vcm) has been computed as (Vcm) = (|xcm(t + t) — 
Xcm(i)|/T)t, with r = 10* and (...) a temporal average. 
O. Biham, A. A. Middleton and D. Levine, Phys. Rev. A 
46, R6124 (1992); T. Nagatani, Phys. Rev. E 48, 3290 

(1993). 

See supplementary material at http: / /link.aps.org/| sup- 

plemental/XXXXX for a movie. 

F.Y. Wu, Rev. Mod. Phys. 54, 235 (1982). 

D.D. Belts, Can. J. Phys. 42, 1564 (1964). 



